home *** CD-ROM | disk | FTP | other *** search
/ MacHack 1998 / MacHack 1998.toast / The Hacks! / PCA Icon Arranger ƒ / Statis.p < prev    next >
Encoding:
Text File  |  1998-06-20  |  7.1 KB  |  315 lines  |  [TEXT/PJMM]
  1.  
  2. unit Statis;
  3.  
  4. interface
  5.     function ChiProb (X: Extended;
  6.                                     F: Integer;
  7.                                     BigX: Boolean;
  8.                                     var Wrong: Boolean): Extended;
  9.  
  10.     function FTest (f: Extended;
  11.                                     df1: Integer;
  12.                                     df2: Integer): Extended;
  13.  
  14.     function Gauss (x: Extended): Extended;
  15.  
  16.     function ggubfs (var dSeed: Extended): Extended;
  17.  
  18.     function FInverse (prob: Extended;
  19.                                     df1: Integer;
  20.                                     df2: Integer): Extended;
  21.  
  22.     function tQuantile (p: Extended;
  23.                                     n: Extended): Extended;
  24.  
  25.     function ZTest (X: Real): Real;
  26.  
  27.  
  28. implementation
  29.  
  30.     function Gauss (X: Extended): Extended;
  31.  
  32.         var
  33.             Y: Extended;
  34.             Z, Z1: Extended;
  35.             W: Extended;
  36.  
  37.     begin { Gauss }
  38.                     (* REFERENCE: CACM            *)
  39.                     (* ALGORITHME: 209            *)
  40.                     (* AUTEUR:D.IBBETSON          *)
  41.                     (* TITRE:GAUSS              *)
  42.                     (* TRADUIT EN PASCAL PAR A.VAUDOR *)
  43.         if X = 0 then
  44.             Z := 0
  45.         else begin
  46.             Y := ABS(X) / 2;
  47.             if Y >= 3 then
  48.                 Z := 1
  49.             else if Y < 1 then begin
  50.                 W := Y * Y;
  51.                 Z := ((((((((0.000124818987 * W - 0.001075204047) * W + 0.005198775019) * W - 0.019198292004) * W + 0.059054035642) * W - 0.151968751364) * W + 0.319152932694) * W - 0.531923007300) * W + 0.797884560593) * Y * 2
  52.             end
  53.             else begin
  54.                 Y := Y - 2;
  55. { Mod. PhC 11/02/98: coupé z en plusieurs morceaux }
  56.                 Z1 := ( ( ( ( ( ( ( ( ( ( ( - 0.000045255659 * Y + 0.000152529290 ) * Y - 0.000019538132 ) * Y - 0.000676904986 ) * Y + 0.001390604284 ) * Y - 0.000794620820 ) * Y - 0.002034254874 ) * Y + 0.006549791214 ) * Y - 0.010557625006 ) * Y + 0.011630447319 ) * Y - 0.009279453341 ) * Y + 0.005353579108 ) * Y - 0.002141268741;
  57.                 Z := ( ( Z1 ) * Y + 0.000535310849 ) * Y + 0.999936657524
  58.             end
  59.         end;
  60.         if X > 0 then
  61.             GAUSS := (Z + 1) / 2
  62.         else
  63.             GAUSS := (1 - Z) / 2;
  64.     end;    (* GAUSS *)
  65.  
  66.  
  67.     function ftest (f: extended;
  68.                                     df1, df2: integer): extended;
  69.         var
  70.             f1, f2, x, ft, vp, xx, theta, sth, cth, sts, cts, a, b, xi, gamma, cbrf: extended;
  71.     begin
  72.                         (* refence:cacm             *)
  73.                         (* algorithme:346             *)
  74.                         (* auteur:j.morris            *)
  75.                         (* titre:ftest probabilities    *)
  76.                         (* traduit en pascal par a.vaudor *)
  77.         if f = 0 then
  78.             ftest := 1.0
  79.         else begin
  80.             f1 := df1;
  81.             f2 := df2;
  82.             ft := 0;
  83.             x := f2 / (f2 + f1 * f);
  84.             vp := f1 + f2 - 2.0;
  85.             if (df1 mod 2 = 0) and (df1 < 500) then begin
  86.                 xx := 1 - x;
  87.                 f1 := f1 - 2;
  88.                 while f1 > 1 do begin
  89.                     vp := vp - 2;
  90.                     ft := xx * vp / f1 * (1 + ft);
  91.                     f1 := f1 - 2;
  92.                 end;
  93.                 if x <> 0 then
  94.                     ft := exp((0.5 * f2) * ln(x)) * (1 + ft)
  95.                 else
  96.                     ft := 0;
  97.             end
  98.             else if (df2 mod 2 = 0) and (df2 < 500) then begin
  99.                 f2 := f2 - 2;
  100.                 while f2 > 1 do begin
  101.                     vp := vp - 2;
  102.                     ft := x * vp / f2 * (1.0 + ft);
  103.                     f2 := f2 - 2;
  104.                 end;
  105.                 if x <> 1 then
  106.                     ft := 1 - exp((0.5 * f1) * ln(1 - x)) * (1 + ft)
  107.                 else
  108.                     ft := 1;
  109.             end
  110.             else if df1 + df2 <= 500 then begin
  111.                 theta := arctan(sqrt(f1 * f / f2));
  112.                 sth := sin(theta);
  113.                 cth := cos(theta);
  114.                 sts := sqr(sth);
  115.                 cts := sqr(cth);
  116.                 a := 0;
  117.                 b := 0;
  118.                 if df2 > 1 then begin
  119.                     f2 := f2 - 2.0;
  120.                     while f2 > 2 do begin
  121.                         a := cts * (f2 - 1.0) / f2 * (1.0 + a);
  122.                         f2 := f2 - 2;
  123.                     end;
  124.                     a := sth * cth * (1.0 + a);
  125.                 end;
  126.                 a := theta + a;
  127.                 if df1 > 1 then begin
  128.                     f1 := f1 - 2;
  129.                     while f1 > 2 do begin
  130.                         vp := vp - 2;
  131.                         b := sts * vp / f1 * (1.0 + b);
  132.                         f1 := f1 - 2;
  133.                     end;
  134.                     gamma := 1.0;
  135.                     f2 := 0.5 * df2;
  136.                     xi := 1;
  137.                     while xi <= f2 do begin
  138.                         gamma := xi * gamma / (xi - 0.5);
  139.                         xi := xi + 1;
  140.                     end;
  141.                     b := gamma * sth * exp(df2 * ln(cth)) * (1.0 + b);
  142.                 end;
  143.                 ft := 1.0 + 0.636619772368 * (b - a);
  144.             end
  145.             else begin
  146.                 f1 := 2.0 / (9.0 * f1);
  147.                 f2 := 2.0 / (9.0 * f2);
  148.                 cbrf := (exp(ln(f) / 3));
  149.                 ft := gauss(-((1.0 - f2) * cbrf + f1 - 1.0) / sqrt(f2 * sqrt(cbrf) + f1));
  150.             end;
  151.             if ft < 0 then
  152.                 ftest := 0
  153.             else
  154.                 ftest := ft;
  155.         end;
  156.     end;   (* ftest *)
  157.  
  158.     function CHIPROB (X: extended;
  159.                                     F: INTEGER;
  160.                                     BIGX: BOOLEAN;
  161.                                     var WRONG: BOOLEAN): extended;
  162.         var
  163.             A, Y, S, E, C, Z: extended;
  164.             EVEN: BOOLEAN;
  165.     begin
  166.         WRONG := FALSE;
  167.         if ((X < 0) or (F < 1)) then begin
  168.             WRONG := TRUE;
  169.             CHIPROB := 0;
  170.         end
  171.         else if X = 0 then
  172.             CHIPROB := 1
  173.         else begin
  174.             A := 0.5 * X;
  175.             EVEN := 2 * (F div 2) = F;
  176.             if ((EVEN) or ((F > 2) and (not BIGX))) then
  177.                 Y := EXP(-A);
  178.             if EVEN then
  179.                 S := Y
  180.             else
  181.                 S := 2 * (GAUSS(-SQRT(X)));
  182.             if F > 2 then begin
  183.                 X := 0.5 * (F - 1);
  184.                 if EVEN then
  185.                     Z := 1
  186.                 else
  187.                     Z := 0.5;
  188.                 if BIGX then begin
  189.                     if EVEN then
  190.                         E := 0
  191.                     else
  192.                         E := 0.572364842925;
  193.                     C := LN(A);
  194.                     while Z < X do begin
  195.                         E := LN(Z) + E;
  196.                         S := EXP(C * Z - A - E) + S;
  197.                         Z := Z + 1;
  198.                     end;
  199.                     CHIPROB := S;
  200.                 end
  201.                 else begin
  202.                     if EVEN then
  203.                         E := 1
  204.                     else
  205.                         E := 0.564189583548 / SQRT(A);
  206.                     C := 0;
  207.                     while Z < X do begin
  208.                         E := E * A / Z;
  209.                         C := C + E;
  210.                         Z := Z + 1;
  211.                     end;
  212.                     CHIPROB := C * Y + S;
  213.                 end
  214.             end
  215.             else
  216.                 CHIPROB := S;
  217.         end;
  218.     end;    (* CHIPROB *)
  219.  
  220.     function ggubfs (var dseed: extended): extended;
  221.         const
  222.             d2p31m = 2147483647.0;
  223.             d2p31 = 2147483648.0;
  224.         var
  225.             x: extended;
  226.             j: integer;
  227.     begin
  228.         x := 16807 * dseed;
  229.         j := trunc(x / d2p31m);
  230.         dseed := x - j * d2p31m;
  231.         ggubfs := dseed / d2p31;
  232.     end;
  233.  
  234.     function FInverse (prob: extended;
  235.                                     df1, df2: integer): extended;
  236.         const
  237.             Epsilon = 10E-6;
  238.         var
  239.             i, j: integer;
  240.             nouveaul, delta, a: extended;
  241.     begin     (* finverse *)
  242.         delta := 1;
  243.         repeat
  244.             delta := delta * 2;
  245.             nouveaul := ftest(delta, df1, df2);
  246.         until nouveaul < prob;
  247.         a := delta;
  248.         while abs(nouveaul - prob) > epsilon do begin
  249.             delta := delta / 2;
  250.             if nouveaul > prob then
  251.                 a := a + delta
  252.             else
  253.                 a := a - delta;
  254.             nouveaul := ftest(a, df1, df2);
  255.         end;
  256.         FInverse := a;
  257.     end;     (* fin finverse *)
  258.  
  259.     function tquantile (p, n: extended): extended;  (* algorithme 396 a.c.m *)
  260.         var
  261.             halfpi, a, b, c, d, x, y: extended;
  262.     begin
  263.         if n = 2 then
  264.             tquantile := sqrt(2 / (p * (2 - p)) - 2)
  265.         else begin
  266.             halfpi := 1.5707963268;
  267.             if n = 1 then begin
  268.                 p := p * halfpi;
  269.                 tquantile := cos(p) / sin(p);
  270.             end
  271.             else begin
  272.                 a := 1 / (n - 0.5);
  273.                 b := 48 / sqr(a);
  274.                 c := ((20700 * a / b - 98) * a - 16) * a + 96.36;
  275.                 d := ((94.5 / (b + c) - 3) / b + 1) * sqrt(a * halfpi) * n;
  276.                 x := d * p;
  277.                 y := exp((2 / n) * ln(x));
  278.                 if y > 0.05 + a then begin
  279.                     x := -1.959963984548082;
  280.                     y := x * x;
  281.                     if n < 5 then
  282.                         c := c + 0.3 * (n - 4.5) * (x + 0.6);
  283.                     c := (((0.05 * d * x - 5) * x - 7) * x - 2) * x + b + c;
  284.                     y := (((((0.4 * y + 6.3) * y + 36) * y + 94.5) / c - y - 3) / b + 1) * x;
  285.                     y := a * sqr(y);
  286.                     if y > 0.002 then
  287.                         y := exp(y) - 1
  288.                     else
  289.                         y := 0.5 * sqr(y) + y;
  290.                 end
  291.                 else
  292.                     y := ((1 / (((n + 6) / (n * y) - 0.089 * d - 0.822) * (n + 2) * 3) + 0.5 / (n + 4)) * y - 1) * (n + 1) / (n + 2) + 1 / y;
  293.                 tquantile := sqrt(n * y);
  294.             end;
  295.         end;
  296.     end; (* tquantile *)
  297.  
  298.     function ZTEST (X: REAL): REAL;
  299.         var
  300.             Y, Z: REAL;
  301.     begin
  302.         X := ABS(X);
  303.         Y := EXP(-X * X / 2) * 0.3989423;
  304.         Z := 1 / (1 + X * 0.2316419);
  305.         Z := 2 * Y * Z * ((((1.330274 * Z - 1.821256) * Z + 1.781478) * Z - 0.3565638) * Z + 0.3193815);
  306.         ZTEST := 1 - Z;
  307.     end;
  308.  
  309. end.
  310.  
  311.  
  312.  
  313.  
  314.  
  315.